Tarea series de tiempo - Carlos Vesga, Santigo Buitrago#
1. Considere la serie de tiempo asociada con los futuros de la criptomoneda Bitcoin desde que comenzó a comercializarse hasta la fecha del día de hoy. Utilice la API de Yahoo Finance para obtener esta serie de tiempo.#
Lectura de datos y EDA#
import warnings
warnings.filterwarnings("ignore")
import yfinance as yf
from yahoo_fin.stock_info import get_data
from datetime import datetime, timedelta
import seaborn as sns
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from numpy import log
import numpy as np, pandas as pd
from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
Warning - Certain functionality
requires requests_html, which is not installed.
Install using:
pip install requests_html
After installation, you may have to restart your Python session.
msft = yf.Ticker("BTC-USD")
msft.info
{'name': 'Bitcoin',
'startDate': 1278979200,
'description': 'Bitcoin (BTC) is a cryptocurrency launched in 2010. Users are able to generate BTC through the process of mining. Bitcoin has a current supply of 19,683,106. The last known price of Bitcoin is 66,797.9989346 USD and is up 3.99 over the last 24 hours. It is currently trading on 11005 active market(s) with $43,228,757,125.25 traded over the last 24 hours. More information can be found at https://bitcoin.org/.',
'maxAge': 86400,
'priceHint': 2,
'previousClose': 63441.96,
'open': 63441.96,
'dayLow': 62850.684,
'dayHigh': 63825.55,
'regularMarketPreviousClose': 63441.96,
'regularMarketOpen': 63441.96,
'regularMarketDayLow': 62850.684,
'regularMarketDayHigh': 63825.55,
'volume': 42898325504,
'regularMarketVolume': 42898325504,
'averageVolume': 34966998847,
'averageVolume10days': 36316462166,
'averageDailyVolume10Day': 36316462166,
'marketCap': 1251204005888,
'fiftyTwoWeekLow': 24797.168,
'fiftyTwoWeekHigh': 73750.07,
'fiftyDayAverage': 66794.59,
'twoHundredDayAverage': 46248.73,
'currency': 'USD',
'fromCurrency': 'BTC',
'toCurrency': 'USD=X',
'lastMarket': 'CoinMarketCap',
'coinMarketCapLink': 'https://coinmarketcap.com/currencies/bitcoin',
'volume24Hr': 42898325504,
'volumeAllCurrencies': 42898325504,
'circulatingSupply': 19683776,
'exchange': 'CCC',
'quoteType': 'CRYPTOCURRENCY',
'symbol': 'BTC-USD',
'underlyingSymbol': 'BTC-USD',
'shortName': 'Bitcoin USD',
'longName': 'Bitcoin USD',
'firstTradeDateEpochUtc': 1410912000,
'timeZoneFullName': 'UTC',
'timeZoneShortName': 'UTC',
'uuid': '74397779-1589-3270-8c45-b7f1a7345b3a',
'messageBoardId': 'finmb_BTC_CCC',
'trailingPegRatio': None}
Podemos ver los datos basicos de la accion, tales como su valor actual, la fecha en que inicio a ser comercializada entre otros.
sns.set_theme()
sns.set_context("paper")
stock = 'BTC-USD'
resolution = '1d'
end_date = datetime.now()
start_date = 1278979200
Como queremos analizar la informacion desde el inicio de la serie de tiempo hasta el momento actual el end_date debe ser igual a “hoy” y el start_date debe coincidir con le fecha de inicio mostrada anteriormente. Como estamos interesados en predecir el valor del día nuestra resolucion es de un día.
def date_format(date_h):
return date_h.strftime('%d/%m/%Y')
BTC_df = get_data(stock, start_date=start_date, end_date=end_date, interval=resolution, index_as_date=False)
BTC_df.head()
| date | open | high | low | close | adjclose | volume | ticker | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2014-09-17 | 465.864014 | 468.174011 | 452.421997 | 457.334015 | 457.334015 | 21056800.0 | BTC-USD |
| 1 | 2014-09-18 | 456.859985 | 456.859985 | 413.104004 | 424.440002 | 424.440002 | 34483200.0 | BTC-USD |
| 2 | 2014-09-19 | 424.102997 | 427.834991 | 384.532013 | 394.795990 | 394.795990 | 37919700.0 | BTC-USD |
| 3 | 2014-09-20 | 394.673004 | 423.295990 | 389.882996 | 408.903992 | 408.903992 | 36863600.0 | BTC-USD |
| 4 | 2014-09-21 | 408.084991 | 412.425995 | 393.181000 | 398.821014 | 398.821014 | 26580100.0 | BTC-USD |
sns.lineplot(data=BTC_df, x=BTC_df.date, y=BTC_df.close);
Al leer y graficar los datos vemos una gran volatilidad en el valor de los BTC, esto es esperado y cabe señalar que entorpece las predicciones realizadas.
import plotly.graph_objects as go
fig = go.Figure(data=[go.Candlestick(x = BTC_df.date,
open = BTC_df.open,
high = BTC_df.high,
low = BTC_df.low,
close = BTC_df.close)
])
fig.update_layout(
title="Bitcoin. (BTC)",
xaxis_title="Day",
yaxis_title="BTC-USD",
font=dict(
family="Courier New, monospace",
size=12,
color="RebeccaPurple"
)
)
fig.update_layout(xaxis_rangeslider_visible=False)
fig
Se grafica tambien un candlestick para facilitar la visualizacion.
Pruebas estadisticas#
Antes de comenzar a pronosticar decidimos realizar algunas pruebas estadisticas para confirmar la viabilidad del modelo ARIMA.
result = adfuller(BTC_df.close)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
ADF Statistic: nan
p-value: nan
Al p-value ser mayor a 0.05 no rechazamos la hipotesis nula, por lo cual la serie estadisticamente no es estacionaria.
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import numpy as np
from statsmodels.tsa.stattools import acf
def plotly_acf_with_confidence(data, lags):
acf_vals, conf_int = acf(data, nlags=lags, alpha=0.05, fft=True) # Obtener ACF y intervalos de confianza
lower_bound = conf_int[:, 0] # límite inferior
upper_bound = conf_int[:, 1] # límite superior
traces = [
go.Scatter(
x=list(range(lags+1)), y=acf_vals, mode='lines', name='ACF',
line=dict(color='blue'),
fill='tozeroy', # Relleno hasta el eje x
fillcolor='rgba(0, 100, 255, 0.2)'
),
go.Scatter(
x=list(range(lags+1)) + list(range(lags+1))[::-1], # x, luego x reversa
y=np.concatenate([upper_bound, lower_bound[::-1]]), # upper, luego lower reversa
fill='toself', fillcolor='rgba(0, 100, 255, 0.3)', line=dict(color='rgba(255,255,255,0)'),
showlegend=False # Ocultar leyenda para este trazo
)
]
return traces
# Suponiendo que tienes 'BTC_df' cargado y listo, y la columna 'close' tiene los precios de cierre
fig = make_subplots(rows=3, cols=2, subplot_titles=("Original Series", "ACF of Original",
"1st Order Differencing", "ACF of 1st Differencing",
"2nd Order Differencing", "ACF of 2nd Differencing"))
# Añadiendo las series y los ACF con intervalos de confianza
fig.add_trace(go.Scatter(x=BTC_df.index, y=BTC_df.close, mode='lines', showlegend=False), row=1, col=1)
fig.add_traces(plotly_acf_with_confidence(BTC_df.close, 3497), rows=[1,1], cols=[2,2])
# Primera diferenciación
first_diff = BTC_df.close.diff().dropna()
fig.add_trace(go.Scatter(x=BTC_df.index[1:], y=first_diff, mode='lines', showlegend=False), row=2, col=1)
fig.add_traces(plotly_acf_with_confidence(first_diff, 3496), rows=[2,2], cols=[2,2])
# Segunda diferenciación
second_diff = BTC_df.close.diff().diff().dropna()
fig.add_trace(go.Scatter(x=BTC_df.index[2:], y=second_diff, mode='lines', showlegend=False), row=3, col=1)
fig.add_traces(plotly_acf_with_confidence(second_diff, 3495), rows=[3,3], cols=[2,2])
# Actualizar layout para mejorar la visualización
fig.update_layout(height=1200, width=1000, title_text="BTC Price Analysis with Interactive ACF Plots")
fig.show()
Graficamos las auto-correlaciones de la grafica original, la de difereciacion de orden 1 y de orden 2. El decaimiento geometrico en la primera grafica nos hace confirmar que la serie de tiempo no es estacionaria, ademas podemos ver que la grafica de diferenciacion de orden 2 entra al lado negativo con rapidez, esto puede ser debido a que esta sobrediferenciada.
2. Repita TODOS los pasos indicados en esta sección para encontrar modelos ARIMA para predecir el precio de Bitcoin con los siguientes horizontes: 7, 14, 21 y 28 días. Utilizar siempre predicciones usando rolling con ventana de predicción continua de un día. Cualquier cantidad de pasos extra para enriquecer su análisis predictivo serán aceptados siempre y cuando sean acordes con lo que indica la teoría de análisis de series de tiempo.#
from statsmodels.tsa.arima.model import ARIMA
import warnings
from statsmodels.graphics.tsaplots import plot_predict
warnings.filterwarnings("ignore")
Hacemos el llamado de la libreria.
def arima_rolling(history, test, orderb):
predictions = list()
for t in range(len(test)):
model = ARIMA(history, order=orderb)
model_fit = model.fit()
output = model_fit.forecast()
yhat = output[0]
predictions.append(yhat)
obs = test[t]
history.append(obs)
return predictions
Definimos la funcion para trabajar en la metodologia rolling forecast.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score
def forecast_accuracy(predictions_list, actual, title_prefix=""):
# Asumiendo que 'actual' es una serie o arreglo que contiene los valores reales
results = []
fig, axes = plt.subplots(nrows=5, ncols=1, figsize=(10, 18))
metrics = ['MAE', 'MSE', 'MAPE', 'RMSE', 'R2']
colors = ['blue', 'green', 'red', 'purple', 'orange']
for idx, forecast in enumerate(predictions_list):
days_predicted = len(forecast)
str_name = f"{title_prefix} {days_predicted} days"
# Calculo de errores
mape = np.mean(np.abs(forecast - actual[:days_predicted]) / np.abs(actual[:days_predicted])) # MAPE
mae = np.mean(np.abs(forecast - actual[:days_predicted])) # MAE
rmse = np.mean((forecast - actual[:days_predicted])**2)**.5 # RMSE
mse = np.mean((forecast - actual[:days_predicted])**2) # MSE
r2 = r2_score(actual[:days_predicted], forecast) # R2
results.append({
'Days': days_predicted,
'MAE': mae,
'MSE': mse,
'MAPE': mape,
'RMSE': rmse,
'R2': r2
})
# Crear un DataFrame para los resultados
df_results = pd.DataFrame(results)
# Graficar cada métrica
for i, metric in enumerate(metrics):
axes[i].plot(df_results['Days'], df_results[metric], label=metric, color=colors[i], marker='o')
axes[i].set_title(f"{metric} by days predicted")
axes[i].set_xlabel('Days')
axes[i].set_ylabel(metric)
axes[i].legend()
plt.tight_layout()
plt.show()
return df_results
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[14], line 4
2 import pandas as pd
3 import matplotlib.pyplot as plt
----> 4 from sklearn.metrics import r2_score
6 def forecast_accuracy(predictions_list, actual, title_prefix=""):
7 # Asumiendo que 'actual' es una serie o arreglo que contiene los valores reales
9 results = []
ModuleNotFoundError: No module named 'sklearn'
Se define esta funcion para el calculo de las metricas de error. Modificamos la funcion enseñada en clase para que acepte listas y grafique cada metrica.
import numpy as np
import scipy.stats as stats
from statsmodels.stats.diagnostic import acorr_ljungbox
def analyze_residuals(predictions_list, actual_data):
results = []
for idx, prediction in enumerate(predictions_list):
# Asumimos que 'actual_data' tiene suficientes datos para comparar con cada 'prediction'
actual_subset = actual_data[:len(prediction)]
# Calcular residuos
residuals = np.array(actual_subset) - np.array(prediction)
# Prueba de normalidad (Shapiro-Wilk)
shapiro_test = stats.shapiro(residuals)
shapiro_pvalue = shapiro_test.pvalue
# Prueba de independencia (Ljung-Box)
# El número de lags recomendado por defecto es min(10, len(residuals) // 5)
lb_lags = min(10, len(residuals) // 5)
lb_test = acorr_ljungbox(residuals, lags=[lb_lags], return_df=True)
# Almacenar resultados
results.append({
'Prediction Length': len(prediction),
'Shapiro P-Value': shapiro_pvalue,
'Ljung-Box P-Value': lb_test['lb_pvalue'].iloc[0]
})
# Imprimir conclusiones
print(f"Residuos para la predicción de {len(prediction)} días:")
print(f"Test de Shapiro-Wilk P-Value: {shapiro_pvalue:.4f}")
if shapiro_pvalue < 0.05:
print("Rechazar la hipótesis nula: Los residuos no son normales.")
else:
print("No rechazar la hipótesis nula: Los residuos parecen ser normales.")
print(f"Test de Ljung-Box P-Value: {lb_test['lb_pvalue'].iloc[0]:.4f}")
if lb_test['lb_pvalue'].iloc[0] < 0.05:
print("Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.")
else:
print("No rechazar la hipótesis nula: Los residuos parecen ser independientes.")
print("\n")
return results
Se crea esta funcion para facilitar las pruebas de normalidad e independencia de los residuales de cada intervalo de prediccion.
n_BTC = len(BTC_df.close); n_test = 28
train_size = n_BTC - n_test
train = BTC_df.close[:train_size]
dates_train = BTC_df.date[:train_size]
test_4w = BTC_df.close[train_size:train_size + n_test]
dates_4w = BTC_df.date[train_size:train_size + n_test]
print("train:", train.shape)
print("test_4w:", test_4w.shape)
train: (3470,)
test_4w: (28,)
Definimos y dividimos el training y el test
train_df = BTC_df[["close"]][:train_size]
test_4w_df = BTC_df[["close"]][train_size:train_size + n_test]
Organizamos el df en train y test dejando los ultimos 28 días para el testing.
best_aic = np.inf
best_order = None
best_mdl = None
pq_rng = range(5)
d_rng = range(3)
for i in pq_rng:
for d in d_rng:
for j in pq_rng:
try:
# print(i, d, j)
tmp_mdl = ARIMA(train, order=(i,d,j)).fit()
tmp_aic = tmp_mdl.aic
if tmp_aic < best_aic:
best_aic = tmp_aic
best_order = (i, d, j)
best_mdl = tmp_mdl
except: continue
Hacemos la hyperparametrizacion del modelo con los parametros moviendose del 1 al 10 usando la metrica aic.
print('aic: {:6.5f} | order: {}'.format(best_aic, best_order))
aic: 56295.88868 | order: (3, 2, 4)
Podemos ver que los mejores parametros segun la metrica aic son (3, 2, 4). Esto es un poco contradictorio con el analisis hecho con las graficas de autocorrelaciones pero trabajaremos con este orden y segun las pruebas estadisticas decidiremos si bajar la diferenciacion a una o mantenerla.
model_aic = ARIMA(train, order=best_order)
model_fit_aic = model_aic.fit()
Entrenamos el modelo con los mejores parametros encontrados.
plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit_aic, 1, ax=ax);
plt.show();
Graficamos el intervalo de confianza del 95% para nuestras predicciones.
test_range = range(4)
predictions_roll_aic = []
for i in test_range:
test = test_4w.head(7*(i+1))
test = test.tolist()
yhat = arima_rolling(train.tolist(), test, best_order)
predictions_roll_aic.append(yhat)
En esta funcion estamos haciendo el rolling forecast para los intervalos de 7,14,21 y 28 días, nos ayudamos en el hecho de ser todos multiplos del 7. Guardamos los resultados de cada pronostico en una lista.
results = analyze_residuals(predictions_roll_aic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.7373
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1296
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.4232
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.3658
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.3687
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5013
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.7656
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5265
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Al utilizar la funcion definida anteriormente podemos ver que todos los supuestos se cumplen para cada residual, por lo que concluimos que el modelo capta la variabilidad de los datos y es posible confiar en las predicciones.
3. Repita el paso 2 ahora sin utilizar rolling. Esto es, realice el pronóstico solo utilizando forecast() para los diferentes horizontes de predicción, 7, 14, 21 y 28 días.#
test_range = range(4)
predictions_for_aic = []
for i in test_range:
# Hacer las predicciones para el siguiente intervalo de tiempo
num_steps = 7 * (i + 1) # Calculando el número de pasos a predecir
yhat = model_fit_aic.forecast(steps=num_steps)
predictions_for_aic.append(yhat)
results = analyze_residuals(predictions_for_aic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.3997
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0792
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.1172
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0402
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.
Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.1362
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1099
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.0445
Rechazar la hipótesis nula: Los residuos no son normales.
Test de Ljung-Box P-Value: 0.0416
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.
Observamos que a diferencia del rolling forecast el supuesto de independencia no se cumple para algunos intervalos de tiempo, se necesita una investigacion en el tema pero para ambitos educativos trabajaremos asi.
4. Realice tablas de error para los ítems 1 y 2, utilizando las métricas: MAPE, MAE, RMSE, MSE, R2. Además, agregue el gráfico de correlación entre la observación real y su predicción en el test.#
Rolling Forecast#
roll_aic = forecast_accuracy(predictions_roll_aic, test_4w, title_prefix="Forecast")
Graficamos y guardamos los resultados conseguidos, por extraño que pareza parece que los resultados mejoran al aumentar el intervalo de prediccion.
Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:
Datos limitados
Sobreajuste
Variacion Aleatoria
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import numpy as np
def plot_prediction_correlation(actual, predictions_list):
# Crear subplots, uno para cada conjunto de predicciones
num_predictions = len(predictions_list)
fig = make_subplots(rows=num_predictions, cols=1,
subplot_titles=[f'Correlacion de {len(pred)} dias predichos' for pred in predictions_list],
vertical_spacing=0.1)
for i, predictions in enumerate(predictions_list):
# Asegurarse de que los datos reales y las predicciones estén alineados
actual_subset = actual[:len(predictions)]
# Calcular coeficientes de correlación
correlation_matrix = np.corrcoef(actual_subset, predictions)
correlation_xy = correlation_matrix[0,1]
r_squared = correlation_xy**2
# Añadir gráfico de dispersión para cada conjunto de predicciones
fig.add_trace(
go.Scatter(
x=actual_subset, y=predictions, mode='markers',
marker=dict(size=8, color='blue', opacity=0.5),
name=f'R^2: {r_squared:.2f}'
),
row=i+1, col=1
)
# Ajustar una línea de regresión lineal
z = np.polyfit(actual_subset, predictions, 1)
p = np.poly1d(z)
fig.add_trace(
go.Scatter(
x=actual_subset, y=p(actual_subset), mode='lines',
line=dict(color='red', dash='dash'),
name='Fit'
),
row=i+1, col=1
)
# Actualizar ejes para cada subplot
fig.update_xaxes(title_text="Valores Reales", row=i+1, col=1)
fig.update_yaxes(title_text="Valores Predichos", row=i+1, col=1)
# Actualizar layout del gráfico
fig.update_layout(
height=300 * num_predictions, width=800,
title_text="Correlacion entre predichos y verdaderos",
showlegend=False
)
fig.show()
Creamos una funcion para graficar las correlaciones entre el valor real y predicho de cada intervalo de tiempo
plot_prediction_correlation(test_4w, predictions_roll_aic)
Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, los mejores R^2 estan para el intervalo de 28 y 14 dias.
Forecast#
for_aic = forecast_accuracy(predictions_for_aic, test_4w, title_prefix="Forecast")
Parece que el mejor intervalo es a 21 días.
Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:
Datos limitados
Sobreajuste
Variacion Aleatoria
plot_prediction_correlation(test_4w, predictions_for_aic)
Todos los intervalos tienen una tendencia positiva,y el mejor R^2 esta a 14 días.
Comparaciones#
import matplotlib.pyplot as plt
def plot_comparative_metrics(df1, df2):
metrics = ['MAE', 'MSE', 'MAPE', 'RMSE', 'R2']
num_metrics = len(metrics)
intervals = df1.index # Asumimos que el índice es el intervalo de días
width = 0.35 # Ancho de las barras
# Crear una figura grande para contener todos los subplots
fig, axes = plt.subplots(nrows=num_metrics, ncols=1, figsize=(10, 6 * num_metrics))
for i, metric in enumerate(metrics):
ax = axes[i] # Acceder al subplot correspondiente
# Graficar las barras para DF1
ax.bar(intervals - width/2, df1[metric], width, label='Rolling', color='blue')
# Graficar las barras para DF2
ax.bar(intervals + width/2, df2[metric], width, label='Forecast', color='orange')
ax.set_title(f'Comparación de {metric}')
ax.set_ylabel(metric)
ax.legend()
axes[-1].set_xlabel('Intervalos de días') # Etiqueta del eje x solo en el último gráfico
axes[-1].set_xticks(intervals) # Asegurar que todos los intervalos estén etiquetados
plt.tight_layout()
plt.show()
plot_comparative_metrics(roll_aic, for_aic)
En terminos generales el rolling forecast es mejor, debido a que minimiza la mayoria de los errores.
5. Repita el análisis desarrollado en los pasos anteriores, considerando ahora el criterio de inferencia Bayesiana (BIC) y el criterio de información de Hannan–Quinn (HQIC) para encontrar el mejor modelo ARIMA y, compare los errores con aquellos obtenidos con el criterio de Akaike.#
BIC#
best_bic = np.inf
best_order = None
best_mdl = None
pq_rng = range(5)
d_rng = range(3)
for i in pq_rng:
for d in d_rng:
for j in pq_rng:
try:
# print(i, d, j)
tmp_mdl = ARIMA(train, order=(i,d,j)).fit()
tmp_bic = tmp_mdl.bic
if tmp_bic < best_bic:
best_bic = tmp_bic
best_order = (i, d, j)
best_mdl = tmp_mdl
except: continue
Volvemos a correr la hyperparametrizacion esta vez optimizando el bic
print('bic: {:6.5f} | order: {}'.format(best_bic, best_order))
bic: 56327.16155 | order: (0, 2, 1)
model_bic = ARIMA(train, order=best_order)
model_fit_bic = model_bic.fit()
plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit_bic, 1, ax=ax);
plt.show();
Graficamos el intervalo de confianza del 95% para nuestro nuevo modelo.
Rolling Forecast#
test_range = range(4)
predictions_roll_bic = []
for i in test_range:
test = test_4w.head(7*(i+1))
test = test.tolist()
yhat = arima_rolling(train.tolist(), test, best_order)
predictions_roll_bic.append(yhat)
Hacemos el rolling forecast para el nuevo modelo
results = analyze_residuals(predictions_roll_bic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.9388
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1408
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.8911
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.3486
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.8613
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5638
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.9871
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.6047
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Cumple todos los supuestos, por lo que concluimos que el modelo capta la variabilidad de los datos y es posible confiar en las predicciones.
roll_bic = forecast_accuracy(predictions_roll_bic, test_4w, title_prefix="Forecast")
Graficamos y guardamos los resultados conseguidos, por extraño que pareza parece que los resultados mejoran al aumentar el intervalo de prediccion.
Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:
Datos limitados
Sobreajuste
Variacion Aleatoria
plot_prediction_correlation(test_4w, predictions_roll_bic)
Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 28 dias.
Forecast#
test_range = range(4)
predictions_for_bic = []
for i in test_range:
# Hacer las predicciones para el siguiente intervalo de tiempo
num_steps = 7 * (i + 1) # Calculando el número de pasos a predecir
yhat = model_fit_bic.forecast(steps=num_steps)
predictions_for_bic.append(yhat)
results = analyze_residuals(predictions_for_bic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.4822
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0864
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.1098
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0435
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.
Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.1307
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1106
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.0391
Rechazar la hipótesis nula: Los residuos no son normales.
Test de Ljung-Box P-Value: 0.0419
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.
Observamos que a diferencia del rolling forecast el supuesto de independencia no se cumple para algunos intervalos de tiempo, se necesita una investigacion en el tema pero para ambitos educativos trabajaremos asi.
for_bic = forecast_accuracy(predictions_for_bic, test_4w, title_prefix="Forecast")
Parece que el mejor intervalo es a 21 días.
Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:
Datos limitados
Sobreajuste
Variacion Aleatoria
plot_prediction_correlation(test_4w, predictions_for_bic)
Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 14 dias.
Comparaciones#
plot_comparative_metrics(roll_bic, for_bic)
En terminos generales el rolling forecast es mejor, debido a que minimiza la mayoria de los errores.
HQIC#
best_hqic = np.inf
best_order = None
best_mdl = None
pq_rng = range(5)
d_rng = range(3)
for i in pq_rng:
for d in d_rng:
for j in pq_rng:
try:
# print(i, d, j)
tmp_mdl = ARIMA(train, order=(i,d,j)).fit()
tmp_hqic = tmp_mdl.hqic
if tmp_hqic < best_hqic:
best_hqic = tmp_hqic
best_order = (i, d, j)
best_mdl = tmp_mdl
except: continue
Hyperparametrizamos optimizando el HQIC
print('HQIC: {:6.5f} | order: {}'.format(best_hqic, best_order))
HQIC: 56313.45958 | order: (3, 2, 4)
model_hqic = ARIMA(train, order=best_order)
model_fit_hqic = model_hqic.fit()
plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit_hqic, 1, ax=ax);
plt.show();
Graficamos el intervalo de confianza del 95%
Rolling Forecast#
test_range = range(4)
predictions_roll_hqic = []
for i in test_range:
test = test_4w.head(7*(i+1))
test = test.tolist()
yhat = arima_rolling(train.tolist(), test, best_order)
predictions_roll_hqic.append(yhat)
Hacemos el rolling forecast para el nuevo modelo
results = analyze_residuals(predictions_roll_hqic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.7373
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1296
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.4232
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.3658
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.3687
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5013
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.7656
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5265
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Cumple todos los supuestos, por lo que concluimos que el modelo capta la variabilidad de los datos y es posible confiar en las predicciones.
roll_hqic = forecast_accuracy(predictions_roll_hqic, test_4w, title_prefix="Forecast")
Graficamos y guardamos los resultados conseguidos, por extraño que pareza parece que los resultados mejoran al aumentar el intervalo de prediccion.
Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:
Datos limitados
Sobreajuste
Variacion Aleatoria
plot_prediction_correlation(test_4w, predictions_roll_hqic)
Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 28 dias.
Forecast#
test_range = range(4)
predictions_for_hqic = []
for i in test_range:
# Hacer las predicciones para el siguiente intervalo de tiempo
num_steps = 7 * (i + 1) # Calculando el número de pasos a predecir
yhat = model_fit_hqic.forecast(steps=num_steps)
predictions_for_hqic.append(yhat)
results = analyze_residuals(predictions_for_hqic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.3997
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0792
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.1172
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0402
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.
Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.1362
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1099
No rechazar la hipótesis nula: Los residuos parecen ser independientes.
Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.0445
Rechazar la hipótesis nula: Los residuos no son normales.
Test de Ljung-Box P-Value: 0.0416
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.
Observamos que a diferencia del rolling forecast el supuesto de independencia no se cumple para algunos intervalos de tiempo, se necesita una investigacion en el tema pero para ambitos educativos trabajaremos asi.
for_hqic = forecast_accuracy(predictions_for_hqic, test_4w, title_prefix="Forecast")
Parece que el mejor intervalo es a 21 días.
plot_prediction_correlation(test_4w, predictions_for_hqic)
Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 14 dias.
Comparaciones#
plot_comparative_metrics(roll_aic, for_aic)
En terminos generales el rolling forecast es mejor, debido a que minimiza la mayoria de los errores.
Comparaciones#
import matplotlib.pyplot as plt
def plot_comparative_metrics_3(df1, df2, df3):
metrics = ['MAE', 'MSE', 'MAPE', 'RMSE', 'R2']
num_metrics = len(metrics)
intervals = df1.index # Asumimos que el índice es el intervalo de días
width = 0.20 # Ancho de las barras, ajustado para acomodar tres modelos
# Crear una figura grande para contener todos los subplots
fig, axes = plt.subplots(nrows=num_metrics, ncols=1, figsize=(10, 6 * num_metrics))
for i, metric in enumerate(metrics):
ax = axes[i] # Acceder al subplot correspondiente
# Graficar las barras para DF1
ax.bar(intervals - width, df1[metric], width, label='AIC', color='blue')
# Graficar las barras para DF2
ax.bar(intervals, df2[metric], width, label='BIC', color='orange')
# Graficar las barras para DF3
ax.bar(intervals + width, df3[metric], width, label='HQIC', color='green')
ax.set_title(f'Comparación de {metric}')
ax.set_ylabel(metric)
ax.legend()
axes[-1].set_xlabel('Intervalos de días') # Etiqueta del eje x solo en el último gráfico
axes[-1].set_xticks(intervals) # Asegurar que todos los intervalos estén etiquetados
plt.tight_layout()
plt.show()
plot_comparative_metrics_3(roll_aic, roll_bic, roll_hqic)
